library(foreign)
library(mediation)
library(ggplot2)
library(plotly)
library(lubridate)
library(reshape2)

sudMed <- read.dta("SudanCEM.dta")
sudMed$`_mi_miss` <- NULL
sudMed$ProvinceName <- factor(sudMed$ProvinceName)

sudImp <- split(sudMed, sudMed$imp)

sudComplete$Education <- NULL
sudComplete$CurrentEcon <- NULL
sudComplete$FutureEcon <- NULL
sudComplete$GovPer <- NULL


sudImp <- sudImp[-1]


#datasets <- list(A1=sudImp$`1`, A2=sudImp$`2`)

datasets <- list(A1=sudImp$`1`, A2=sudImp$`2`, A3=sudImp$`3`, A4=sudImp$`4`, A5=sudImp$`5`, A6=sudImp$`6`,
                 A7=sudImp$`7`,A8=sudImp$`8`, A9=sudImp$`9`, A10=sudImp$`10`, A11=sudImp$`11`, A12=sudImp$`12`, 
                 A13=sudImp$`13`, A14=sudImp$`14`, A15=sudImp$`15`, A16=sudImp$`16`, A17=sudImp$`17`, 
                 A18=sudImp$`18`, A19=sudImp$`19`, A20=sudImp$`20`, A21=sudImp$`21`, A22=sudImp$`22`, 
                 A23=sudImp$`23`, A24=sudImp$`24`, A25=sudImp$`25`, A26=sudImp$`26`, A27=sudImp$`27`, 
                 A28=sudImp$`28`, A29=sudImp$`29`, A30=sudImp$`30`, A31=sudImp$`31`, A32=sudImp$`32`,
                 A33=sudImp$`33`, A34=sudImp$`34`, A35=sudImp$`35`, A36=sudImp$`36`, A37=sudImp$`37`,
                 A38=sudImp$`38`, A39=sudImp$`39`, A40=sudImp$`40`)

datasetsSens <- list(sudImp$`1`, sudImp$`2`, sudImp$`3`, sudImp$`4`, sudImp$`5`, sudImp$`6`,
                     sudImp$`7`,sudImp$`8`, sudImp$`9`, sudImp$`10`, sudImp$`11`, sudImp$`12`, 
                     sudImp$`13`, sudImp$`14`, sudImp$`15`, sudImp$`16`, sudImp$`17`, 
                     sudImp$`18`, sudImp$`19`, sudImp$`20`, sudImp$`21`, sudImp$`22`, 
                     sudImp$`23`, sudImp$`24`, sudImp$`25`, sudImp$`26`, sudImp$`27`, 
                     sudImp$`28`, sudImp$`29`, sudImp$`30`, sudImp$`31`, sudImp$`32`,
                     sudImp$`33`, sudImp$`34`, sudImp$`35`, sudImp$`36`, sudImp$`37`,
                     sudImp$`38`, sudImp$`39`, sudImp$`40`)

mediators <- c("Freeriding")

outcome <- c("CritGov2")

treatments <- rep("ArabSpring", 40)

#treatments <- c("ArabSpring", "ArabSpring")

covariates <- c("Rural + Female + lnincome + Age + ProvinceName")

weights <- c("cem_weights")

test <- mediations(datasets, treatments, mediators, outcome, covariates, families=c("gaussian","gaussian"), weights, interaction=FALSE, conf.level=.95, sims=1000)

output <- amelidiate(test)

testSens <- mediate(output$model.m, output$model.y, sims=1000, treat = "ArabSpring", mediator = "Freeriding")
testSens.c <- medsens(testSens, rho.by = 0.05)

summary(output)
plot(output)


#sensitivity
med.fit <- lm(Freeriding ~ ArabSpring + Rural + Female + lnincome + Age + ProvinceName, weights = cem_weights, data=sudImp$`3`)
out.fit <- lm(CritGov2 ~ Freeriding + ArabSpring + Rural + Female + lnincome + Age + ProvinceName, weights = cem_weights, data=sudImp$`2`)
med.out <- mediate(med.fit, out.fit, treat="ArabSpring", mediator="Freeriding", robustSE=TRUE, sims=1000)
sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims=1000)
summary(sens.out)
plot(sens.out, sens.par = "rho", main = "Criticize Government", ylim = c(-0.2, 0.2))

imputedSens <- list()
imputedSens <<- append(imputedSens, list(sens.out))
imputedSens[1]
plot(imputedSens[1], sens.par = "rho", main = "CritGov2", ylim = c(-0.2, 0.2))

imputed.mediation <- function(x){
  imputedSens <<- list()
  rhoSens <<- c()
  for (i in x){
    med.fit <- lm(Freeriding ~ ArabSpring + Rural + Female + lnincome + Age + ProvinceName, weights = cem_weights, data = i)
    out.fit <- lm(CritGov2 ~ Freeriding + ArabSpring + Rural + Female + lnincome + Age + ProvinceName, weights = cem_weights, data = i)
    med.out <- mediate(med.fit, out.fit, treat="ArabSpring", mediator="Freeriding", robustSE=TRUE, sims=1000)
    sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims=1000)
    rhoAtZero <- sens.out$err.cr.d
    imputedSens <<- append(imputedSens, list(sens.out))
    rhoSens <<- append(rhoSens, rhoAtZero)
  }
}

imputed.mediation(datasetsSens)
rhoAnalysis <- data.frame(rhoSens,listOne=rep(1, 40))
mean(rhoSens)
library(ggplot2)
ggplot(data = rhoAnalysis, aes(x = rhoSens, y=listOne)) +
  geom_bar(stat = "identity") +
  xlab(bquote(rho)) + ylab("Count") + labs(caption = bquote(bar(rho) == -0.35)) +
  ggtitle(bquote(rho ~ "at ACME = 0 for Imputed Datasets"))